Scaling theory for spontaneous imbibition in random networks of elongated pores 



Zcinab Sadjadi 1 and Hciko Rieger 1 

'Theoretical Physics, Saarland University, 66041 Saarbriicken (Germany) 
(Dated: October 18, 2012) 

We present a scaling theory for the long time behavior of spontaneous imbibition in porous media 
consisting of interconnected pores with a large length-to- width ratio. At pore junctions the meniscus 
propagation in one or more branches can come to a halt when the Laplace pressure of the meniscus 
exceeds the hydrostatic pressure within the junction. We derive the scaling relations for the emerging 
arrest time distribution and show that the average front width is proportional to the height, yielding 
a roughness exponent of exactly ft = 1/2 and explaining recent experimental results for nano-porous 
Vycor glass (NVG). Extensive simulations of a pore network model confirm these predictions. 
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The dynamics of imbibition front of an invading fluid 
in disordered media has attracted substantial scientific 
attention, from statistical physics [U-IH to material sci- 
ence [Hi. Besides its scientific interest, understanding the 
mechanisms of imbibition in a porous matrix is of impor- 
tance in industrial processes such as oil recovery, food 
processing, impregnation, chromatography, and agricul- 
ture SSQ. 

During imbibition the liquid-gas interface advances 
and broadens. The time evolution of the invading front 
follows simple scaling laws, which are independent of the 
micro-structure and the details of the fluid (341^. rem- 
iniscent of the universality of critical phenomena. Var- 
ious physical aspects are involved in the imbibition of 
a liquid inside a porous matrix, such as viscous drag, 
capillarity, gravity and volume conservation. The often 
complex topology of the porous matrix induces local fluc- 
tuations in capillary pressures at the interface as well 
as hydraulic permeabilities in the bulk. Despite these 
complexities, the average position of the front H during 
a purely spontaneous imbibition evolves as if(t)^ 1 / 2 , 
known as Lucas- Washburn law 0, [3, EH . This scaling 
behavior is valid down to nanoscopic pore scales EMI- 

While the invading front exhibits a common slow- 
broadening dynamics for a wide range of materials 0- 
lllj ]. the results of recent experiments on nano-porous 
Vycor glass (NVG) reveal that the roughening dynamics 
might depend on the micro-structure [19j. The elonga- 
tion of pores, quantitatively described by their length- 
to-width radius, appears to play an important role and 
two extreme limits can be distinguished: (i) short pores 
with comparable length and diameter. In materials like 
paper, sand, randomly packed glass beads, etc., where 
the pore space is highly interconnected j8l-tLU|. neigh- 
boring menisci coalesce, a continuous imbibition front 
forms and an effective surface tension emerges. Due to 
the latter menisci advancement is spatially highly corre- 
lated [20[, which reduces the height fluctuations of the 
front by limiting menisci advancement beyond the aver- 
age front position and drawing forward the menisci lag- 
ging behind. This forms a continuous liquid-gas interface 



and smoothens the front, (ii) elongated pores. Other 
porous materials like rock, soil, and porous glasses con- 
sist of sponge-like topologies with reduced connectivity 
and elongated pores [5, 2M, 22 1- For example, NVG is a 
silica substrate with an interconnected network of long 
cylindrical pores with characteristic radii of 3—5 nm. In 
Ref. an anomalously fast interface roughening has 
been observed, representing a new universality class for 
spontaneous imbibition, emerging for large pore aspect 
ratio. Here, the interface is not able to establish an effec- 
tive surface tension, leading to strong height fluctuations 
of the menisci. 

In this letter we present a scaling theory for sponta- 
neous imbibition in porous media consisting of a network 
of interconnected elongated pores (Fig JT]) . It is based on 
the observation that at pore junctions the meniscus prop- 
agation in the branch with the larger radius can come to 
a halt when the Laplace pressure of the meniscus exceeds 
the hydrostatic pressure within the junction. This leads 
to the emergence of voids behind the invasion front and 
concomitantly to anomalously fast front broadening as 
observed experimentally in NVG [l9[ . It is predicted that 
the distribution of the meniscus arrest times scales with 
the square of the height of the meniscus, which implies 
that the ratio of the average invasion front width and 




FIG. 1: Sketch of a junction (a) in a pore network with elon- 
gated pores (b). r* and Pl,% denote the radius and Laplace 
pressure, respectively, in pore i, and po denotes the hydro- 
static pressure in the junction. In (b) (h(t)) and (Pl) denote 
the average height at time t and the average Laplace pressure. 
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the average front height is independent of time. This 
implies that roughening is maximal with an exponent 
j3 = 1/2, establishing a universality class different from 
those known before for spontaneous imbibition Q. We 
then test these predictions in extensive simulations of a 
pore network model. 

We analyze spontaneous imbibition of a wetting liq- 
uid in a porous medium similar to porous glasses, which 
consists of a network of elongated pores with a length- 
to- width ratio of the order of 10, i.e. elongated, cylinder- 
like pores with random radii interconnected at pore junc- 
tions as sketched in Fig. 1. The bottom pores are con- 
nected to a liquid reservoir with pressure p — 0. We 
assume that in each pore a liquid-gas interface forms, de- 
noted as meniscus, that gives rise to a Laplace-pressure 
PL =—2cr/r, where a is the surface tension of the liquid 
and r the pore radius. If the pore radii vary between r m i n 
and r max , the average radius is denoted by (r). Then, 
on large scales, the average height is expected to vary 
as d/dt{h(t)) = — (Pi) I (h(t)) , which implies the Lucas- 
Washburn law (h(t)} cx t 1/2 . 

Consider now a junction at height ho, where a pore 
branches into two (see Fig. la). One branch has radius 
r\, the other r 2 > r%, yielding the Laplace pressures 
Pl,i= —2a/ri. Let Pq be the hydrostatic pressure within 
the junction. As long as Pl,i > Po the meniscus in 
branch 2 is arrested. In the following we will answer 
the question how long the meniscus in branch 2 will be 
arrested and we will implicitly assume that it does not 
get annihilated by the filling of the pore from its other 
end. This means that we assume the radius r 2 also to be 
larger than the radius of the other branch of the junction 
of the other end. This reduces only the probability of 
this event by a ^-dependent factor. 

-Po = Po(t) is a function of time and depends on how 
far the front has propagated and can be estimated as fol- 
lows: Let the average front height be (h(t)). On average 
one expects the bulk pressure to decrease linearly from 
bottom to top: 

P((h(t)))/Po = (h(t))/h (1) 

Therefore, with P((h(t))) = (P L ) = -2a(l/r) the av- 
erage Laplace pressure, one obtains Po = — 2<r(l/r) • 
ho/(h(t)} and the condition P = Pl,2 f° r the arrested 
meniscus to resume propagation (at time Resume) reads 

(MWnnc)) = h Q r 2 {l/r) . (2) 

This equation has far reaching consequences: 
1) The larger r 2 the longer the meniscus is arrested, and 
the average height that the front has to reach before 
the meniscus resumes propagation is proportional to the 
height where it stopped with a proportionality constant 
larger than one. 2) The time r for which the meniscus is 
arrested is proportional to the time t a t op , when it stopped 

t oc t stop . (3) 



To see this we note that with ([2]) one has (h(t stop +r)) = 
h(t s t op )r 2 (l/r). With Lucas- Washburn (h(t stop + r)) oc 
(tstop +t) 1/2 and assuming that h(t stop ) oc ^/ 2 p , too, for 
the relation between the height and the time when the 
considered meniscus stopped, one obtains ([3]). 

3) Consequently from ([3]) 

r oc h 2 (t st0 p ) = hi;, (4) 

which implies that the probability distribution of arrest 
times for menisci arrested at height h will scale as 

Ph (r) = h- 2 p(r/h 2 ) . (5) 

4) The height difference w (t lcsumo ) = (h(t lcsumc )) - ho 
is a measure for the local width of the propagation front 
(at the lateral coordinates of the position of the arrested 
meniscus) at time t esumc . The ratio of this local width 
and the average height is ^ trca """j = 1 — (^(l/r)) -1 , 
which is independent of the time t resume . Thus all ar- 
rested menisci will contribute a time independent amount 
to the ratio of the average width w(i) and average height. 
Since the width cannot grow faster than h(t) this implies 

w(t)/(h(t)) = const. , (6) 

implying wit) oc t 1 / 2 , i.e. a roughening exponent @ = 1/2 
Note that the invasion front dynamics is now expected 
to be be completely determined by the meniscus arrests, 
which in turn depend exclusively on the pore radii dis- 
tribution and the height dependent hydrostatic pressure. 
Consequently one expects no lateral correlations in the 
meniscus heights to emerge, as observed in fl9| . 

The scaling theory presented here neglects all geomet- 
ric and topological details of a pore network. To test 
its predictions, in particular the strongest ([5J and ([5]), 
we analyzed the following microscopic model for sponta- 
neous imbibition in a pore network with elongated pores 
[lil I23I f24"j : A two-dimensional square lattice of cylindri- 
cal capillaries inclined at 45° is considered, which con- 
sists of N x and N y nodes in horizontal and vertical direc- 
tions, respectively. Capillaries, interconnected at nodes, 
have the same length L and random radii uniformly dis- 
tributed over [r av — 6, r av + S] . The average aspect ratio 
2r av /i is set to 5. The pressure at the bottom nodes at- 
tached to the liquid reservoir is set to zero, the pressure at 
a moving meniscus is the Laplace pressure. Here we ne- 
glect gravity, which is justified as long as capillary forces 
are much larger than gravitational forces 2a /r ^> pN y L, 
where p the specific weight of the liquid. This is the case 
for instance in experiments with NVG [25| . 

The hydrostatic pressures at the nodes of the network 
drives the dynamical evolution of the menisci configu- 
rations. To calculate the temporal change of the fill- 
ing heights in the partially filled capillaries, one needs to 
know the node pressures which themselves depend on the 
menisci configuration. The node pressures Pi are deter- 
mined by the boundary conditions plus the conservation 
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FIG. 2: The mechanisms of menisci advancement in the 
pore-network model (a) after reaching a node, and (b) during 
backward motion and arrest of a meniscus due to negative 
pressure difference. 



of volume flux at each node: Y]j Qj—O, which is equiv- 
alent to Kirchhoff's law. Here, Q\ is the volume flux 
flowing from node i into the capillary j attached to it. 
The sum runs over all of the four capillaries of node i 
and is valid for all wet nodes in the system. According 
to Hagcn-Poiseuille's Law Q\= - 4 AP? /h{, with 



d-. 



=7r(r^) 4 /8?7 and AP?=Pi-P 3 Li . Here, 



are the radius, the length and the Laplace pressure of 
the meniscus in capillary j of node i, respectively, and rj 
is the viscosity of the liquid. By numerically solving the 
resulting set of linear equations we compute Pi and thus 
Qj . These arc then inserted into the equation of motion 
for the heights given by Q^=7r(r^) 2 dh\/dt. To integrate 
these differential equations an implicit Euler scheme with 
variable time step At is employed giving the new posi- 
tions h\ . When a meniscus reaches the end of a capillary 
it immediately moves an infinitesimal distance Sc^O.OlL 
into the adjacent capillaries, creating new menisci, as 
shown in Fig. [Ha). This avoids the microscopic treat- 
ment of the filling of the junction and is valid as long 
as the filling time of the node is negligible. The filling 
time was estimated in [27j and is indeed orders of magni- 
tudes smaller than the meniscus arrest times as long as 
capillary forces are much larger than gravitational forces. 
When two menisci meet, they vanish, thus the capillary is 
entirely filled. If, due to a negative pressure difference, a 
meniscus retracts, it proceeds backward as long as its dis- 
tance from the back node is larger than 6 [see Fig. HJb)]. 
When it reaches <5, the meniscus is stuck there until the 
driving pressure difference is again positive. During this 
arrest time, the pressure calculation is modified with the 
corresponding capillary being blocked. We made sure 
that the simulation results we present in the following 
are independent of the choice of S. 

First we examined the arrest events during the front 
propagation. Fig.rSJa) shows three snapshots of the prop- 
agating and arrested menisci in the invasion front at three 
different times. As can be seen with increasing height 
more and more menisci are arrested. By counting the 
number of menisci arrested at height h for a time r we ob- 
tained the arrest time distribution ph{r), which is shown 
in Fig. EJb). It broadens systematically with increasing 
height h and the inset shows that it scales nicely with h 2 
as predicted by eq.© ph{t) = h~ 2 p(rjh 2 ). This seal- 
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FIG. 3: (a) Snapshots of the arrested (blue circles) and ad- 
vancing (red circles) menisci in the invasion front at three dif- 
ferent times. Broken (full) lines represent empty (full) pores, 
H is the average height at the corresponding time, (b) The 
probability distribution ph(r) of the arrest times r of menisci 
arrested at height h. Inset: Scaling plot according to eq.([5|. 
(c) Average arrest time (r) at a given height h, dashed line 
proportional to h 2 , see eq.Q. Inset: Probability for a menis- 
cus to be arrested at height h. 5 r /r = 0.1, N x = 8, 100 
disorder realization. 



ing relation implies that the average arrest time (r) at 
height h should be proportional to h 2 ( cq.2]). FigurelSJc) 
shows (r) as a function of h in a log-log plot in compari- 
son with the straight line for h 2 . For large heights small 
but systematic downward deviations occur which might 
be due to the fact that ph(j) appears to have an alge- 
braic tail (see[3fb)), implying that rare events (large r) 
would dominate the average over disorder realizations. 
This might then be underestimated for a restricted num- 
ber of samples. Finally the inset of Fig. (3jc) shows that 
the probability for a meniscus to be arrested quickly ap- 
proaches one for increasing h: At sufficiently long times 
the invasion front mainly consists of arrested menisci. 

Next we compute the average height H(t)= (h(t)) and 
width w(t)=({h 2 (t)) - (hit)) 2 ) 1 / 2 of the imbibition front . 
Figure 0] shows that the ratio w(t)/H(t) approaches a 
constant value for large times, which confirms the pre- 
diction dH]). Together with H(t) oc A 2 this implies that 
asymtotically the width also increases as w(t) oc t 1 / 2 , 
implying a roughness exponent (3 = 1/2. The initial de- 
crease of w(t)/H(t) indicates a pre-asymptotic increase 
of w(t) with an exponent slightly smaller than 1/2, as re- 
ported in [l||. The data shown in Fig. |4] are independent 
of N x confirming the absence of finites size effects and 
lateral correlations in the meniscus heights. 
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porous medium consisting of random networks of elon- 
gated pores. We tested the predictions in extensive sim- 
ulations of a pore network model. Meniscus arrest times 
at pore junctions are shown to scale with the age of the 
invasion front whose width is therefore proportional to its 
average height. This establishes a universality class for 
invasion front broadening that is realized in nano-porous 
Vycor glass [3] and is expected to determine roughen- 
ing dynamics in similar porous media. Since meniscus 
arrest is solely determined by the relation of radii of the 
pores emanating from one junction, it should be possi- 
ble to relate dynamical quantities accessible via light or 
Neutron scattering to characteristics of the pore radius 
distribution of the porous medium. 



FIG. 4: (a): The width of imbibition front w(t) (open sym- 
bols) and the ratio of width to height w(t)/H(t) (full symbols) 
vs. time for different disorder on radius. The lattice size is 
N x xN v = 16xl000. (b): The same plot for different lateral 
size of the system, and <5t.=0.1. The dashed line indicates the 
slope /3=l/2 as a guide to the eye. 



The invasion front thus involves a finite fraction of 
the occupied volume and comprises connected clusters 
of empty pores whose size distribution gets broader with 
increasing time. Based on the conditions for meniscus ar- 
rests presented above one can derive a scaling form for the 
distribution of cluster sizes as follows. Consider an empty 
cluster that contains S pores at time t (i.e. with width 
w(t)). Its lateral size scales as C ~ S x / df and its surface 
area as J- ~ S d "' df , where df and d$ are the bulk and 
surface fractal dimension of the empty clusters (df = d 
and d s = d — 1 in the case of compact clusters). Almost 
all pores in the boundary T of the cluster have arrested 
menisci, and for a meniscus to be arrested the radius of 
its pore has to be larger than the radius of an adjacent 
pore, which is an event that occurs with some probability 
q < 1. Assuming that the conditions for meniscus arrest 
in all boundary pores are independent from one another 

pores, 



Acknowledgements: We thank P. Huber and D.-S. 
Lee for stimulating discussions. 



d s /d f 



and the boundary consists of the order of S' 
the probability for collective meniscus arrests in bound- 
ary pores is proportional to cxp(— a • S ds ^ df ), where the 
constant a involves Inq and a geometric factor. Since 
the lateral dimension S 1 ^ 1 of the cluster must not ex- 
ceed the width w(t), one obtains for the probability of 
an arbitrary empty pore in the front region to belong to 
a connected cluster with S empty pores: 

q s = A/" -1 S exp(-a • S^' d f) • ~g(S 1 ' d i /w(t)) , (7) 

where J\f is a normalization factor and g(x) is a scaling 
function that is 1 for i <C 1 and for x — > 1. A cluster 
analysis of our simulation of the 2c? pore network model 
confirms the stretched exponential behavior of qs at large 
times (w(i) > S 1 /^) with d s /d f close to 0.5 = (d-l)/d. 

In conclusion we have presented a scaling theory for 
the imbibition of an arbitrary wetting liquid through any 
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